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Abstract 



Results of a phenomenological Monte Carlo calculation for a 2D electron- 
phonon Holstein model near a Peierls-CDW transition are presented. Here 
the zero Matsubara frequency part of the phonon action is dominant and we 
approximated it by a phenomenological form that has an Ising-like Peierls- 
CDW transition. The resulting model is studied on a 32x32 lattice. The 
single particle spectral weight A{k,uj)^ the density of states N{u;), and the 
real part of the conductivity cri{uj) all show evidence of a pseudogap which 
develops in the low-energy electronic degrees of freedom as the Peierls-CDW 
transition is approached. 
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I. INTRODUCTION 



In the electron-phonon Holstein model, local phonon modes with frequency c^o are cou- 
pled to the electron charge density on each site of a lattice. Here we consider this model 
on a 2D square lattice with a near-neighbor, one-electron transfer matrix element t. At 
half-filling this system undergoes an Ising-like phase transition to a Peierls-CDW phase at a 
transition temperature Tc. We are interested in developing a more detailed picture of what 
happens to the electronic degrees of freedom as this phase transition is approached. 

When the temperature is reduced, the lattice site displacement-displacement correlation 
length increases, diverging as (T/Tc — 1)^^ on an infinite lattice. This behavior clearly alters 
the usual Migdal, electron-phonon renormalized, quasi-particle description [|T|. In particular, 
as Tc is approached, a pseudogap in what were the quasi-particle degrees of freedom develops. 
Previous studies of this type of model have used various analytic approximations as 
well as quantum Monte Carlo methods @J3l- However, the non-gaussian nature of the 
lattice displacement fluctuations and their strong coupling to the electronic degrees of 
freedom have raised questions about the analytic approximations and the relatively small, 
2D lattices (10 x 10) that have been studied with exact Monte Carlo algorithms have limited 
their range. Here we investigate this problem using a static, phenomenological Monte Carlo 
approach, that allows us to treat the non-gaussian fluctuations and strong coupling effects 
on larger lattices (32x32) in the pseudogap regime. In a static approximation, only the 
zero frequency Matsubara part of the effective phonon action is kept, which corresponds 
to treating the thermal lattice fluctuations classically. These classical thermal fluctuations 
determine the critical behavior at the Peierls-CDW transition and the long-wave-length, 
low-frequency response of the lattice. 

Previously, Monte Carlo simulations of classical Heisenberg antiferromagnetic spins in- 
teracting with electrons have been used to model the mangenites and recently a classical 
xy action has been used to study the effect of phase fluctuations on a d-wave BCS system 
TO| , p]l|] . Monien has studied the cross-over from Gaussian to non-Gaussian behavior for a ID 



electronic system coupled to classical incommensurate CDW fluctuations p|. A static field 
approximation in which the action was obtained by integrating out the fermions has also 
been used to treat the positive U Hubbard model |]T2[. Here we will approximate the effective 
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static part of the phonon action using a phenomenological Ising form. We are interested in 
the effect that an Ising-hke Peierls-CDW transition has on the electronic properties. We find 
that the single particle spectral weight A{k,u), the density of states N{u), the spin suscep- 
tibility, and the real part of the conductivity (Ji{uj) all show evidence of a pseudogap which 
develops in the low-energy electronic degrees of freedom as the Peierls-CDW transition is 
approached. 

The Holstein Hamiltonian, which we will study, consists of a one-electron hopping term 
t, local Einstein phonons of frequency ujq and an interaction which couples the on-site local 
phonon displacement Xj and the site charge density. 

H = -tT. (4c,o + clc,^) + E [Tb + ^ ^i)+9j2 (1-1) 

{ij)<T i \ ) ia 

The sum in the first term is over the set of nearest neighbor sites on a 2D square lattice. 
The c|g. fermion operator creates an electron of spin a on site i. Pi and Xj are the conjugate 
momentum and position operators for the local phonon mode on the z*'^ site, g is the electron- 
phonon coupling, and = c\^Ciu. 

In the usual quantum Monte Carlo approach, the fermions are traced over leaving an 
expression for the partition function in terms of an effective action ^[{^^(r)}] for the lattice 
displacement field Xj(r). Here we imagine that this has been done and then express Xj(r) 
in terms of its Matsubara frequency transform. 

a;,(r)=^x,(a;Je-^--^ (1.2) 

m 

with {ujm = 2m7[T). As the transition temperature is approached, the electron-phonon 
coupling leads to a softening of the phonon frequencies and for Q = (7r,7r), ujq vanishes 
at Tc. In the region near where ujq <^ 27rT, the effective action is dominated by the 
zero frequency Matsubara part. As discused by Bickers and one of the authors |]12| for the 
Hubbard model, one could now proceed to directly calculate the = contribution to the 
action by integrating out the fermions in the presence of a given static field configuration. 
Here, however, we will approximate the effective action by an Ising-like form 

s[{x.}] = ( ^ E + c E ^^^A (1-3) 
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with Xi = Xi{uJm = 0). Here, a, 6, and c vary slowly with temperature and we will treat them 
as phenomenological constants. In the following we will set the scale of Xi so that a = 1 and 
take 6 = 2 and c = 1/8 in units where the hopping t = 1. These values imply {xf) ~ 1 and 
give a mean field transition temperature T^'^^ = 0.5t. 

Now the calculation of the electronic properties reduces to the problem of electrons 
moving in a static random external field {xi} with 

H = -t J2 {4aCj<^ + h-c) + 9 ^XiTli^ . (1.4) 
{i,j)cr ia 

The probability distribution of the external field is proportional to 

P[{xi}] = ^"P^'t^^"-^]^ (1.5) 



with S'[{xi}] given by eq. (|T]3|) and 



Z = j {^dx)j eM-S[{xi\]) . (1.6) 

In the following we will set the coupling g in eq. ( |1.4| ) equal to 1. 

Even on a 32 by 32 lattice, certain physical quantities such as the density of states show 
significant size effects due to the degeneracies of the simple nearest-neighbor tight-binding 
spectrum we are using. In order to break these degeneracies and obtain results closer to the 



continuum limit, we follow Assaad [|13| and introduce a small magnetic field that vanishes 
in the thermodynamic limit. The magnetic field is introduced via the Peierls phase factors 
and the electronic Hamiltonian becomes 



H = {^JijCja + h.c}j + g Yl^i^i^y = cLH^j[{x}] c^„ (1.7) 

(i.j)u ia- (ij)a- 



with 



t 



-te *o (1.8) 



where B = V x A and $0 is the flux quantum. The size of the magnetic field is B = ^q/N 
where N = L x L[L = 32) is the number of lattice sites. We use a symmetric gauge with 



m = ^{-iyj.,0) (1.9) 



and boundary conditions such that 



As Assaad showed, the introduction of such a magnetic field removes the level degeneracy 
spreading the single particle states over the bandwidth. This procedure provides a significant 
reduction in the finite size effects and we will use it for all the quantities that will be 
calculated. While it breaks translational invariance for a finite lattice size L, the invariance 
is naturally restored when L goes to infinity. On our finite lattice, we will simply take the 
average over all lattice sites of the quantities of interest. 



II. CALCULATIONS 

We are interested in the fermion spectral function and density of states which can be 
obtained from the single particle Green's function: 

G,,.(r,rO = -{{Tr{cUr)cUr')})f){.} (2.1) 

The double average ((■ ■ is defined as 

_ Tr[exp(-/?g) ■ ■ ■] 
^ " Tr[exp{-/3H)] ^ " ^ 

= {Ud^r)exp{-S[{x}])... (2.3) 



where H, S[{x}] and Z are given by eqs. ( |1.4|) , (|1.3| ), ([L6| ) respectively. In order to carry 
out the fermion trace it is convenient to make the canonical transformation {cj^, Ci^} — > 
{ipicryi^ia} that diagonalizes the Hamiltonian matrix eq. ( |1.4[ ). 

c.. = EO».[{x}]^.. ; cl = J20l[{x}]^l (2.4) 

r r 

j:OU{^mA{x}]0,r[{x}] = e,.[{x}]6rs (2.5) 

ij 

Y.0U{x}]0„.[{x}] = 6rs (2.6) 

After carrying out the fermion trace, and analytically continuing to real frequencies, the 
retarded single particle Green's function is obtained as 

gU^) = (j:Orr[{x}]o*A{x}] Arn^)i . (2.7) 

u — er[{x\\ + lO+ ' {x} 



While averaging over the Monte Carlo configurations {x} restores the translational symmetry 
breaking associated with the individual displacement field configuration, the introduction 
of the magnetic field, eq. (|1.8| ), breaks translational symmetry. As noted, we are willing to 
accomodate this feature in order to reduce the finite size effects that appear in the frequency 
dependence of A{k, uj) and N{uj). With this in mind, we define the Fourier transform of the 
single particle Green's function as: 

G"^(k,^) = ^Ee^'^'''"''^^Gg.H (2-8) 
where N is the number of sites. The spectral function A{]s.,uj) is then given by: 

TT 

^^^.ic(R,-R,) ^0,.[{a;}]0;j{x}]5(.. - ^r[{x}]))^^^ (2.9) 

ij r 

and the single-particle density of states N{uj) is 

TT iV ^ 



j^T.m{x}]\'6{u;-er[{x}])\ (2.10) 

/ {x} 



In the numerical calculations, the 5-functions are then replaced by Lorentians of width 
r] = OMt: 

- e.) - (2.11) 

TT [UJ^ - try + T 

We are also interested in the real part of the conductivity cri(ci;) which can be obtained 
from the analytic continuation of the usual current- current correlation function 

A,(r) = ((T^Jf (r)j;(0)});)|.} . (2.12) 

Here, 

Ji ~ {ti+xiCi^/j-ijCifj tii^xCia-Ci+xfrj (2.13) 

with X a unit displacement in the x-direction. Analytically continuing the Matsubara fre- 
quency transform of Dij{T) to real frequencies, we obtain the retarded propagator whose 
imaginary part is given by 



with 

^ii:rs[{a;}] = mi-r sm*.rs + '^j;sr"^*.,^ - rUi-rsmj-sr " m*.r^ml,^ (2.15) 

and 

m^jr-s = tii+£Oi+£r[{x}]0*,[{x}] (2.16) 

The optical conductivity o"i(co') is then obtained from the imaginary part of the retarded 
current-current correlation function ImD^(ci;), eq. ( p.l4|) as 



^iM = ^ElmA» (2.17) 

A check on the calculation of (Ji{uj) is provided by the f-sum rule relating the average 
kinetic energy per site in the x-direction to the integral over all frequencies of cri(co'). The 
kinetic energy operator for site i is 



Its average is readily obtained as 

(<^^)/)w = (ERe(t..+A[{a:}]0,;^,,,[{a;}]/(6.))\ (2.19) 
leading to the f-sum rule 

poo YJ- 

duja,{uj) = --«K^>f>{xy (2.20) 
In the same manner, one can study the magnetic spin susceptibility 

X(g,r) = ((T{M,(g,r)M,(g,0)});). (2.21) 

with 

MM-Z--'-'C-^^) (2-22) 
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Here we are interested in the q goes to zero static Pauli spin susceptibility which is given by 
Xo = - (I E 0.r[{x}]OU{x}]0,,[{x}]0*A{x}] J^^'J^^ Z ) (2.23) 



{x} 



The calculations were carried out on a 32x32 lattice. The probability distribution P[{x}], 
eq. ( |1.5D , is sampled with a Hybrid Monte Carlo algorithm |T^-|IB|, which we have found 
to be more efficient than the local Metropolis updating scheme. A step of the Hybrid 
Monte Carlo algorithm consists of a proposed move and an accept/reject procedure. A 
move is proposed by generating Gaussian distributed auxiliary momenta with zero mean 
and variance equal to the temperature. Newton's equations of motion using S[{xi}]/ [3, eq. 
( |1.3D , as the potential energy are integrated for n steps with a time-reversible integration 
scheme (necessary for detailed balance), such as the Verlet algorithm that we used. The 
proposed move is accepted with probability mm(l, exp(— /^Ai?)) where /\E is the difference 
in energy between the final and initial configuration and the energy is the sum of the potential 
energy S'[{xj}]//9 and the kinetic energy of the auxiliary momenta. The accept/reject step 
corrects for the approximate integration of Newton's equations of motion and yields an 
exact algorithm. The number of molecular dynamics steps n and the time step used for 
the integration must be chosen to maximize the efficiency of the sampling procedure. We 
set these parameters by maximimizing the average distance between an initial configuration 
and the configuration after one step of the Hybrid Monte Carlo procedure. We found that 
20 molecular dynamics steps with a time step of 0.1 gave a consistently good performance 
for all of the temperatures of interest. For these values of the parameters, the acceptance 
ratio, which depends on how well the energy is conserved by the approximate integration of 
Newton's equations, turns out to be nearly independent of temperature and between 85% 
and 90% for the 32 x 32 lattice. Note that this acceptance ratio does depend on the size 
of the system, because the energy is an extensive quantity, and decreases as the system size 
increases. 

Other parameters of the algorithm that must be set are the number of thermalization 
steps taken before any measurements and how many steps of the algorithm are necessary to 
generate statistically independent measurements. While these parameters can be different 
for different observables, we have looked at the average value of the action, eq. ( |1.3| ), and 
the order parameter Xg = ^J2i{~^yxi which we expect to be the slowest observable to 
thermalize and decorrelate. Since for a finite system one can sample configurations with 
both Xg > and < during a run, one must actually monitor the absolute value of Xs 
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in the simulation. Looking at the Monte Carlo time series of values of \xs\ and S'[{a;i}], we 
found that 10,000 thermalization steps of the Hybrid Monte Carlo algorithm were adequate 
to equilibrate the system for our 32 x 32 lattice. 

The number of steps of the sampling procedure necessary to generate statistically in- 
dependent measurements was determined by the standard data blocking procedure. The 
values of the observables are averaged over 1,2,3, ...n consecutive Hybrid Monte Carlo steps 
and the variance computed. When the number of steps n is greater than the correlation time 
of the sampling procedure, the variance no longer changes. A plot of the variance versus 
n allows one to determine the number of steps needed to generate statistically independent 
measurements easily. For all temperatures of interest, 100 steps of the Hybrid Monte Carlo 
algorithm were adequate to decorrelate the values of \xs\ and S'[{xj}]. 

Note that the parameters of the algorithm need not be the same for each temperature. 
The correlation and thermalization times are smaller at high temperatures. But given that 
by far the most time consuming part of the simulations is the calculation of the electronic 
observables, any fine tuning of the Hybrid Monte Carlo sampling procedure would only lead 
to marginally small gains in computer time. 

We carried out simulations at 11 temperatures starting with 0.5t progressively reducing 
the temperature down to 0.25t. Except for the highest temperature, the initial configu- 
ration was taken as an equilibrium configuration of the previous simulation. Each of the 
runs consisted of 1024 sample configurations with 100 Hybrid MC steps between successive 
measurements. The calculations were carried out on 64 processors of an IBM SP parallel 
computer, using a different seed of the random number generator on each processor. Thus 
16 sample configurations from each processor were needed and the total length of a run 
on a given node was 1600 Hybrid Monte Carlo steps. Close to and below Tc, this is not 
sufficiently long for the system to sample configurations with positive and negative values of 
Xs during a run and thus the system is trapped near one of the two peaks {xs > or < 0) 
of the probability distribution P[{x}], eq. (|1.5|). Since one expects the system to end up in a 
state with Xs > or < with equal probability as the temperature is lowered, on average 
one ought to sample configurations with Xs > on half of the processors and configurations 
with Xs < on the other half. 
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The data analysis was carried out with the multiple histogram technique of Ferrenberg 
and Swendsen [|T^. This technique uses the information available from all of the runs and 
allows one to obtain values for the observables for any temperature in the range covered by 
the simulations. But it requires the range of energies sampled from simulations at neighbor- 
ing temperatures to overlap and this determined the number of runs (11) that we carried 
out. The statistical error estimates were obtained with the bootstrap method. Our simu- 
lation data was resampled 128 times, which allows an estimation of the error bars to a bit 
better that 10%. The multiple histogram calculation of the partition functions and observ- 
ables must be repeated for each bootstrap sample and the statistical errors obtained as the 
variance of the 128 multiple histogram values of the observables. With 1024 samples, the 
statistical errors were typically of the order of 1% in the spectral function A{k, u) and 0.3% 
in the density of states N{uj) and the optical conductivity <Ji{uj). 



III. RESULTS 

For the parameters a = 1,6 = 2, and c = 1/8 that we have chosen, the mean field 
transition temperature T^^^ = 0.5t. Results for the Ising-specific heat on a 32x32 lattice 
are plotted in Fig. 1 and imply that the actual transition temperature is close to = 
0.27t. As discussed, we are interested in the temperature region T*^^ > T > and 
particularly the behavior as T approaches where the uj^ = Matsubara frequency part 



of the action is dominant. In this region, the Peierls-CDW rms gap A^.^^ = {x^) is 
of order t for the parameters we have selected. Thus, the basic coherence length, which 
varies as t/Arms is of order the lattice spacing and the important length scales are the 
Peierls-CDW correlation length ^pcdw{T) and the quasi-particle thermal correlation length 
C,th{T). For a non-interacting system, the thermal quasi-particle length for k ~ (7r/2,7r/2) 



is ^thiT) ^ 2t/7rT, while for k ~ (vr,0), ^th(T) varies as 2yt/T/7r due to the Van Hove 
singularity. 

Results for C,th{T) for k = (7r/2,7r/2) and k = (tt, 0) are shown as the long- and short- 
dashed lines respectively in Fig. 2. The temperature dependence of the Peierls-CDW corre- 
lation length C,pcDw{T) obtained from a Monte Carlo simulation of (xiXj) using the action 



given by eq. ( |1.3|) is also plotted in Fig. 2 [|T8|. Here, one sees that as the temperature is 
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lowered from T^^ towards Tc, the Peierls-CDW correlation length exceeds ^th(T) as Tc is 
approached. We also see that this occurs at a higher temperature for quasi-particles moving 
in the [1,0] direction (kp = (vr,0)) than for quasi-particles moving in the [1,1] direction 
(kp = (7r/2,7r/2)) as expected. 

Turning now to the electronic properties, we first look at the single particle density of 
states N{u!) given by eq. ( p.lO|) . Results for N{u!,T) versus u for various values of the 
temperature are shown in Fig. 3 and a three-dimensional plot of N{uj,T) versus u and T 
is shown in Fig. 4. The lattice fluctuations at the higher temperatures have wiped out the 
Van Hove peak. As the temperature is decreased and the Peierls-CDW correlation length 
exceeds the quasi-particle thermal correlation length, coherence peaks appear in N{uj,T). 
This behavior is similar to what has been found for the 2D Hubbard model near half-filling 
as T goes towards zero and the antiferromagnetic correlation length exceeds the thermal 
quasi-particle length [|l^,|3 ■ 

In Figs 5a and 5b we show the single-particle spectral weights A{k,uj,T) for k = 
(7r/2,7r/2) and (vr, 0) respectively. Again, at higher temperatures, the thermal lattice fluc- 
tuations give rise to a broad peak in these spectral weights, due to the static random nature 
of the lattice displacements. As the temperature is lowered and the lattice displacement 
correlation length increases, coherent quasi-particle peaks, associated with the Peierls-CDW 
phase which onsets at T^. are clearly seen. A careful examination shows that these coherence 
peaks appear at a higher temperature for k = (vr, 0) than for k = (7r/2,7r/2). Thus, the 
Peierls-CDW gap appears to open first at the k = (vr, 0) region of the Fermi surface and 
then at a lower temperature for the diagonal k = (7r/2,7r/2) region. This simply reflects 
the fact that the Peierls-CDW correlation length first exceeds the k = {tt, 0) quasi-particle 
thermal correlation length, as seen in Fig. 2. 

In Fig. 7, we plot the static spin susceptibility Xo(T) versus T. As the temperature 
decreases towards Tc, Xo(T) clearly shows the development of a pseudogap. Finally, turning 
to the real part of the conductivity, we have plotted ai{uj) versus u and T in Fig. 7. As 
noted, the area under ai^uj) should be equal to — 7r/2 times the x- kinetic energy per site. If 
the integral of cri(co') over u is performed analytically for each Monte Carlo configuration, 
this sum rule is satisfied to within statistical errors. The thermally randomized lattice 
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gives rise to elastic scattering at higher temperatures. As the temperature is lowered, the 
pseudogap in the quasi-particle spectrum gives rise to a shift in ai{uj) spectral weight from 
low frequencies to frequencies above 2 Arms- 



IV. CONCLUSION 

The electronic properties of a 2D electron-phonon model which undergoes a commensu- 
rate Peierls-CDW transition have been studied within a static approximation which takes 
into account only the ujm = Matsubara fluctuations of the lattice. Using a Monte Carlo 
calculation based upon a phenomenological Ising-like form for the lattice displacements, the 
effect of these classical thermal lattice fluctuations on the one-electron density of states, 
the single particle spectral weight, the spin susceptibility and the optical conductivity were 
determined. Each of these quantities shows evidence of a depletion of low-energy electronic 
degrees of freedom pseudogap opens. 

For the parameters which we have chosen, ,^0 ~ t/^rms is of order the lattice spacing 
so that the important lengths for this 2D lattice are the thermal quasi-particle length and 



the Peierls-CDW correlation length [^,^. As the temperature is lowered toward T^, the 
Peierls-CDW correlation length flrst exceeds the thermal quasi-particle correlation length 
for electrons near k ~ (tt, 0) and a pseudogap begins to open over the parts of the fermi 
surface associated with the Van Hove singularity of the non-interacting system. Then as the 
temperature is further decreased, a temperature is reached such that ^pcdw(^) exceeds the 
thermal quasi-particle coherence length for the diagonal {tt/2,tt/2) region and a pseudogap 
is opened over the entire fermi surface. 

This continuous depletion of the single particle density of states at the fermi energy 
appears as a suppression of the Pauli spin susceptibility. In addition, spectral weight in 
the optical conductivity is shifted from lower frequencies to higher frequencies, eventually 
appearing as a peak above Aq, resembling cri(ci;) in the ordered Peierls-CDW phase. 

While the results we have presented are in agreement with various theoretical and nu- 
merical expectations, we believe that this numerical method which takes into account the 
uJm = Matsubara fluctuations represents a useful approach to a more detailed understand- 
ing of the electronic properties of a 2D electron-phonon system as a commensurate second 
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order phase transition is approached. Prom a diagramatic point of view, it sums all Feyn- 
man graphs associated with cUm — phonon propagators. The finite Matsubara frequency 
phonons renormalize the parameters a, b, and c in the action, but it is the cUm = fluctu- 
ations that dominate the behavior as T ^ T^. By treating a, b, and c phenomenologically, 
we are able to run the simulation on substantially larger lattices than previous quantum 
Monte Carlo simulations of the Holstein model. This has allowed us a more detailed look 
at the manner in which the electronic properties are effected as the Peierls-CDW transition 
is approached. 
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FIGURES 
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FIG. 1. Monte Carlo results for the Ising-specific heat versus T for a 32 x 32 lattice are shown 
for a = 1, 6 = 2, and c = 1/8. All of our results are for these same parameters. For this figure, we 
used 11 runs of 100,000 measurements each. The statistical errors are of the order of the width of 
the line. 



15 



6 




I I 

0.25 0.3 0.35 0.4 0.45 0.5 

T/t 

FIG. 2. The Peierls-CDW correlation length ^pcDw{T) versus T is shown as the solid curve. 
The non-interacting thermal, quasi-particle, one-electron correlation lengths for k = (7r/2,7r/2) 
and k = (tt, 0) are shown as the long and short dashed curves respectively. 



16 



0.3 



0.2 



0.1 







T = 0.25t 
T = 0.30t 
T = 0.35t 
T = 0.40t 
T = 0.45t 
T = 0.50t 




2 3 4 5 



-5 -4 -3 -2-10 1 

(o/t 

FIG. 3. The single-particle density of states N{uj) versus u> is shown for various temperatures. 



the statistical errors (not shown) are of the order of the linewidth 
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FIG. 4. A plot of A'^(<^, T) versus lo and T. The statistical errors which aren't shown for clarity 
are of the order of 0.3% 
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Spectral Function at k = (7r/2,7r/2) 
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Spectral Function at k = (7t,0) 
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FIG. 5. a) The single-particle spectral weight A{k,uj) versus uj and T for k = (7r/2,7r/2) and 
b) for k = (tt,0). The statistical errors are of the order of 1% 
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FIG. 6. The static spin susceptibility xo(9 = 0, a; = 0) versus T showing the opening of a gap 
as the temperature is lowered. 
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Optical Conductivity 




FIG. 7. Plot of the real part of the optical conductivity ai{uj,T) versus lo and T showing the 
redistribution of spectral weight to higher frequencies as Tc is approached. The statistical errors, 
not shown, are of the order of 0.3% 
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